Damage Detection as Inverse Problems for Distributed

system (2.1) can thus be written as h  w(t); i+ hA2(q) _ w(t); i+ hA1(q)w(t); i = hf(t; q); i w(0) = w0 ; _ w(0) = w1 (2.2) for all 2 V . Adopting somewhat commonly employed abbreviated notation, we let h ; i denote h ; iV ;V throughout the remainder of this discussion. The admissible parameter set Q is assumed to lie in a metric space ~ Q with metric d while the operators Ai(q) are assumed to satisfy the following regularity hypotheses uniformly in q 2 Q: (1) Symmetry: hA1(q) ; i = hA1(q) ; i; (2) Boundedness: jhAi(q) ; ij cij jV j jV with ci independent of q; (3) Coercivity: hAi(q) ; i kij jV with ki independent of q; (4) Continuity with respect to parameters: jh(Ai(q) Ai(~ q)) ; ij d(q; ~ q)j jV j jV . Under these assumptions, one can establish well-posedness of (2.1) (or equivalently, (2.2)) and then use it as the de ning dynamical system in inverse or parameter estimation problems. (In fact, one can establish well-posedness under much weaker assumptions{see [10]). One seeks to estimate the parameters q (e.g., coe cients in some DPS written abstractly as (2.1)) from dynamic observations of the system (2.1). This raises the important consideration as to what will be measured in the dynamic experiments from which we obtain our observations. In mechanical experiments, there are a number of popular measurement devices including accelerometers (yielding acceleration at a point on the structure), laser vibrometers (velocity), proximity probes (displacements), strain gauges (strain), and piezoceramic patches (accumulated strain) in smart structures. In thermal experiments one can use infrared sensors to provide measurements of temperature. For all such measurement devices, the resulting observations can be employed in a general least squares output formulation of the parameter estimation problem. In such cases, the problems are stated in terms of nding parameters which give the best t of the parameter-dependent solutions of the partial di erential equations to dynamic system response data collected with various excitations. The general least squares parameter estimation problem can be formulated as follows. For a given discrete set of measured observations z = fzigNt i=1 corresponding to model observations zob(ti) at times ti as obtained in most practical cases, one considers the problem of minimizing over q 2 Q the least squares output functional J(q; z) = ~ C2 n ~ C1fw(ti; ; q)g fzigo 2 ; (2.3) where fw(ti; ; q)g are the parameter dependent solutions of (2.1) evaluated at each time ti; i = 1; 2; : : : ; Nt and j j is an appropriately chosen Euclidean norm. Here the operators ~ C1 and ~ C2 are observation operators that depend on the type of observed or measured data that is available. The operator ~ C1 may have several forms depending on the type of sensors being used. For example, when the collected data zi consists of time domain displacement, velocity, 3 or acceleration values at a point y, the operator ~ C1 involves di erentiation (either 0, 1, or 2 times, respectively) with respect to time followed by pointwise evaluation at ti and y. While the most commonly encountered least squares problems involve time domain data, it is often advantageous to t the data in a frequency domain setting. To treat these alternate possibilities, one may introduce a second observation operator ~ C2. This operator ~ C2 may be the identity (corresponding to time domain identi cation procedures as already described) or may be related to the Fourier transform (corresponding to identi cation in the frequency domain). If the identi cation is carried out in the frequency domain and the operator ~ C2 is a Fourier transform related operator, then an appropriate cost functional is Ĵ(q; z) = Nf X̀=10@ 1jfkẁ(q) fkz̀ j2 + 2 N X̀ j= n`fjW (kẁ + j; q)j jZ(kz̀ + j)jg21A : (2.4) HereW (k; q) and Z(k) are the Fourier series coe cients of ~ C1fw(ti; y; q)g and fzig respectively, and fkẁ and fkz̀ are the (kẁ)th vibration frequency of the solution W (k; q) and the (kz̀)th frequency of the observation data Z(k). Moreover, 1; 2 are weighting constants, and n`; N` are certain lower and upper limits associated with the width (or the support) of the `th spike. In formulating (2.4), one assumes that there are a nite and distinct number Nf of nontrivial \spikes", i.e., vibration frequencies or \signi cant modes" among the Z(k) and the number of nontrivial spikes of the solution W (k; q) is the same as Nf . (Further details may be found in [10, 12].) The minimization in these general abstract parameter estimation problems involves an in nite dimensional state space and possibly an in nite dimensional admissible parameter set (of functions). To obtain computationally tractable methods, one must consider some type of approximations in the context of the variational formulation (2.2). Let V N be a sequence of nite dimensional subspaces of V , and QM be a sequence of nite dimensional sets approximating the parameter set Q. We denote by PN the orthogonal projections of V onto V N . Then a family of approximating estimation problems with nite dimensional state spaces and parameter sets can be formulated by seeking q 2 QM which minimizes JN (q; z) = ~ C2f ~ C1fwN (ti; ; q)g fzigg 2 ; (2.5) where wN (t; q) 2 V N is the solution to the nite dimensional approximation of (2.2) given by h  wN (t); i+ hA2(q) _ wN(t); i+ hA1(q)wN(t); i = hf(t; q); i wN (0) = PNw0 ; _ wN (0) = PNw1; for 2 V N . Under reasonable (and veri able) assumptions on Q; QM and V N (e.g., compactness of Q; QM and convergence of QM to Q; V N to V as N;M ! 1 in an appropriate sense see [3, 7, 10] for detailed statements), one can show that solutions q and qN;M to the minimization problems for (2.3) and (2.5), respectively, exist. Moreover, qN;M ! q in the metric of Q as N;M !1. 4 3 Thermal NDE Methods One application of the methodology outlined in the previous section is to the problem of detection and characterization of damages using thermal probes ([4, 5, 6]). To illustrate we consider a thin 2-D domain G(q) as depicted in Figure 1. In this geometry the sides and 0 are assumed known while the back side @G(q) is the damaged area that is not observable directly. y x 0 1 Γ0 Γ Γ ` @G(q) G(q) Figure 1. Thin 2-D sample with back surface damage In experiments, the front surface 0 is heated and a time record of the temperature on this surface is obtained using an infrared sensor. The sides are assumed to be perfectly insulated while the damaged section is assumed to disrupt heat ow in a manner consistent with an insulated surface. Thus heat ow in the sample is described by the following model in strong form ( is the thermal di usivity, is the thermal conductivity): @u @t 4 u = 0 on G(q) u(0) = u0 on G(q) @u @n = 0 on [ @G(q) @u @n = f on 0: The associated inverse problems consist of estimating @G(q) from observations of temperature u on 0. Also, the input f on 0 is typically not directly measurable and must also be estimated from the surface temperatures. Both the surface @G(q) and the input ux are assumed to be parameterized by some vector parameter q and we attempt to estimate a best value of q so that the model provides solutions with values on 0 close to those observed in the experimental data. 5 A least squares t of the time dependence of the temperature of the front surface of the sample can be used to attempt to determine the back surface geometry. This requires solving a partial di erential equation (PDE) for shapes calculated from the parameters of the t. Since no analytical solution is known for this case, a nite element technique can be used to solve the PDE. One approach to this problem would be to vary the nite element grid as a function of the parameters describing the shape of the objects. This approach however is computationally di cult and may lead to discontinuities in the derivatives of the solutions with respect to the parameters (an unpleasant situation for most optimization codes). A second approach involves using the method of mappings which is computationally much simpler. The method of mappings transforms the problem of solving a relatively simple PDE on a complex geometry to a coordinate system where the complexity is shifted to the PDE and the geometry is relatively simple. Details of the application of this technique to this type of problem are given in [6]. In the example described here, the shape of the sample was mapped into a coordinate system where the shape of the sample was rectangular. The sample here is taken to be a thin plate with a rectangular geometry, with the exception of a portion of the back surface. The variation in the back surface shape is assumed the same for the thickness of the plate; therefore the heat ow is two dimensional. The shape of the back surface is assumed to be given by r(x) with a constraint on the shape requiring that r(0) and r(1) are both equal to `. The case of r(x) equal to ` for all x corresponds to a rectangular geometry for the plate (i.e., no damage). We assume that the initial temperature u0 of the sample is known. The variational form of this problem can be written as ZG(q) @u @t + ru(t) r ! dV = Z 0 f ds : Here G(q) is the two dimensional plate volume, 0 is the boundary corresponding to y = 0, denotes the trace operator on 0 and is any member of the class of test functions used in formulating the variational equation (see [6] for details). This problem can be transformed to a coordinate system where the shape of the plate was rectangular using the transformation x = x2, and y = y2r(x2)=`. For this coordinate system the variational form is given by the expression (again see [6]) ZV2 " @~ u @x2 @ ~ @x2 r0y2 r @~ u @x2 @ ~ @y2 + @~ u @y2 @ ~ @x2!+ (r0)2y2 2 + `2 r2 @~ u @y2 @ ~ @y2# dV2 + ZV2 "@~ u @t r0 r @~ u @x2 + (r0)2 r2 y2 @~ u @y2# ~ dV2 = Z 0 f ~ r̀ds where ~ u and ~ refer to the temperature and test functions in the transformed coordinate system, respectively, and V2 is the rectangular region 0 x2 1; 0 y2 `. As can be seen by a comparison of equations, the transformation of coordinate system transfers the complexity of the problem to the PDE. However, computationally it is easier to solve the transformed equation than to solve a problem involving variation of the grid to account for variation in shape of the domain. 6 To carry out the least-squares optimization in this problem, the curve r(x) and the input f were parameterized by a nite dimensional vector parameter (again call it q). In this case, linear interpolation was used and best t piecewise linear r; f were obtained. A rst order in time version of the methodology outlined in Section 2 was used to provide theory and computational methods. Bilinear splines were used for the Galerkin approximations of the thermal dynamics. Complete details including the successful identi cation of damage from experimental data are given in [6]. 4 Numerical Methods for Beam Damage Detection To illustrate the damage detection methodology of Section 2 in the context of dynamic structural models with variable coe cients, we consider here a damaged cantilever beam. For this discussion, it is assumed that a single pair of piezoceramic patches with edges at y1; y2 are mounted to the surface of a beam of length ` as depicted in Figure 2. The thickness of the beam and patches are denoted by h and hpe, respectively, while both beam and patches are assumed to have width b. Furthermore, it is assumed that the beam is damaged by a circular hole of radius r and center yc. In this section, numerical techniques satisfying the hypotheses of Section 2 will be discussed to illustrate the issues which must be considered when implementing PDE-based damage detection methods. A Galerkin-based spline approximate will be compared with a standard Hermite nite element method to emphasize the issues. While the development is speci c to the beam example, similar issues must be addressed when developing numerical methods for implementing PDE-based damage detection techniques for complex structures. c w(t,y) y y y 1 2 Figure 2. Cantilever beam with surface-mounted piezoceramic patches and a hole of radius r centered at yc. 7 4.1 Strong Form of Model As detailed in [2, 10], moment and force balancing (Newton's laws) can be used to obtain a strong form of the modeling equations ~ (y)  w(t; y) + _ w(t; y) + hf cD(y) _ w00(t; y) +g EI(y)w00(t; y)i00 = f(t; y) (4.1) with the essential boundary conditions w(t; 0) = w0(t; 0) = 0 (4.2) at the xed end and the natural boundary conditions (f cD(y) _ w00(t; y) +g EI(y)w00(t; y))jy=` = 0 (f cD(y) _ w00(t; y) +g EI(y)w00(t; y))0jy=` = 0 at the free end. In this Euler-Bernoulli model, w denotes the transverse beam displacement, f denotes distributed surface forces on the beam and _ w;w0 denote temporal and spatial derivatives, respectively. Due to the presence of the patches and damage, the structural density, sti ness and Kelvin-Voigt parameters ~ (y);g EI(y) and g cDI(y) are spatially varying and have the forms ~ (y) = hb+ 2 pehpeb pe(y) SA(y) d(y) g EI(y) = E 1 12h3b+ Epe 2 3a3b pe(y) ESI(y) d(y) g cDI(y) = cD 1 12h3b+ cDpe23a3b pe(y) cDSI(y) d(y) : The patch and damage components are localized to their respective regions by the characteristic functions pe(y) = ( 1 y 2 [y1; y2] 0 otherwise ; d(y) = ( 1 y 2 [yc r; yc + r] 0 otherwise : Furthermore, a3 (h=2 + hpe)3 (h=2)3 while E;Epe, and cD; cDpe are the Young's moduli and Kelvin-Voigt damping coe cients of the beam and the piezoceramic patches, respectively. The general shape functions SA and SI indicate the area and moment of inertia of the missing region. For the speci c case of a circular hole with center yc and radius r, these shape functions can be represented as SA(y) = hS(y), SI(y) = 1 12h3bS(y) where S(y) = 2qr2 (y yc)2 : The reader is referred to [2, 10] for discussion of more general shape functions. Furthermore, if the force to the beam is provided solely through input of a voltage u(t) to the patches, f is given by KB 00 pe(y)u(t) where KB = (h+ hpe)bEped31 (see [8]). Due to the discontinuities in the physical parameters and control inputs which are due to the patches, classical solutions (w(t; y) 2 H4(0; `)) cannot be obtained to (4.6). To alleviate such di culties and to reduce smoothness requirements on basis functions, it is advantageous to consider a weak form of the modeling equations derived via Hamilton's principle. 8 4.2 Weak Form Via Hamilton's Principle We illustrate rst the conservative case (no damping or force). With T;U denoting the kinetic and potential energy, respectively, Hamilton's principle states that dynamics of the beam yield a stationary value of the action integral A[w] = Z t1 t0 (T U)dt ; for arbitrary t0; t1 when compared with admissible variations in the motion [17]. This requirement can be expressed as d d"A[w+ " ]j"=0 = 0 : (4.3) For the beam under consideration, the kinetic and potential energies are given by T (t) = 1 2 Z ` 0 ~ (y)( _ w(t; y))2 dy U(t) = 1 2 Z ` 0 g EI(y)(w00(t; y))2 dy : The admissible variations are assumed to be of the form = (t) (x) where (t1) = (t2) = 0 and 2 V = H2 L(0; `) fv 2 H2(0; `) : v(0) = v0(0) = 0g. Expansion of the functional in (4.8) and integration by parts in time then yields the weak formulation Z ` 0 ~ (y)  w +g EI(y)w00 00 dy = 0 which must hold for all 2 H2 L(0; `). Inclusion of damping contributions and external forces and/or moments can be accomplished via an extended Hamilton's principle. For the damped beam subject to forces generated through piezoceramic patch inputs, this yields the weak form Z ` 0 h~ (y)  w + _ w +g EI(y)w00 00 + f cD(y) _ w00 00i dy = Z y2 y1 KB 00dy u(t) (4.4) for all 2 V . It should be noted that smoothness requirements on the solutions w are reduced and this makes the weak form appropriate for consideration of the beam (as well as more general structures) with embedded or surface-mounted actuators. For cases in which parameters and control inputs are su ciently smooth to yield w 2 H4(0; `), integration by parts can be used to demonstrate the equivalence of the strong form (4.6) and weak form (4.9). To place the model in the framework of Section 2, we take H = L2(0; `) and V = H2 L(0; `) as de ned above. With the inner products h ; iH = Z ` 0 ~ dy h ; iV = Z ` 0 g EI 00 00 dy ; 9 the embedding V ,! H is dense and continuous. The duality product h ; iV ;V is then the extension by continuity of the inner product h ; iH from V H to V H. The operators Ai(q) : V ! V , i = 1; 2, are de ned by (A1(q) )( ) = Z ` 0 g EI 00 00 dy (A2(q) )( ) = Z ` 0 g cDI 00 00 dy : With the input operator B 2 L(U; V ), where U is the input space, given by hBu(t); iV ;V = Z y2 y1 KB 00 dy u(t) ; the system (4.9) is equivalent to the operator equation (2.1) in V or the weak form (2.2) with f(t) = Bu(t). The operators A1; A2 de ned in this manner satisfy the symmetry, boundedness, coercivity and parameter continuity conditions described in Section 2. When combined with the approximation techniques discussed next, the results of Section 2 yield a theoretically rigorous and computationally tractable method for estimating the unknown structural parameters which include damage criteria. 4.3 Galerkin Spline Approximation One technique for approximating the solutions to (4.9) is through cubic B-spline expansions in a Galerkin setting. The basis in this case is constructed from canonical cubic B-splines f~ ig of the form illustrated in Figure 3 (see [16, pages 78-80] for details regarding the de nition of these splines). Because the canonical splines satisfy the relations ~ i(xi) = 4, ~ i(xi 1) = ~ i(xi+1) = 1 and ~ 0i(xi) = 0; ~ 0i(xi 1) = 1; ~ 0i(xi+1) = 1 when mapped to the interval [xi 2; xi+2], the choice 1(y) = ~ 0(y) 2~ 1(y) 2~ 1(y) i(y) = ~ i(y) ; i = 2; ; N + 1 yields a basis f igN+1 i=1 which satis es the essential boundary conditions (4.7). One then de nes the approximating space V N =spanf ig V and constructs approximate solutions through the expansion wN (t; y) = N+1 X j=1 wN j (t) j(y) : (4.5) This approximation technique satis es the necessary convergence criteria alluded to in Section 2 and is computationally e cient and accurate. It should be noted that through construction, the splines have continuous rst and second derivatives. By incorporating all smoothness constraints in the basis, the problem is reduced from that of nding a solution which satis es mixed interpolatory and smoothness constraints to one which is purely interpolatory. This will be in contrast to nite element techniques in which degrees of freedom arise from mixed constraints. 10 −2 −1 0 1 2 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 3. Canonical cubic spline ~ i(x) on the reference interval [ 2; 2]. Orthogonalization of the residual with respect to basis functions then yields the approximating ODE system MN  WN (t) + CN _ WN (t) +KNWN (t) = FNu(t) where the mass, sti ness and forcing components are de ned by MN = 2664 R0̀ e 1 1dy R0̀ e N+1 1dy ... ... R0̀ e 1 N+1dy R0̀ e N+1 N+1dy 3775 ; KN = 2664 R0̀ g EI 00 1 00 1dy R0̀ g EI 00 1 00 N+1dy ... ... R0̀ g EI 00 N+1 00 1dy R0̀ g EI 00 N+1 00 N+1dy 3775 FN = 2664 R y2 y1 KB 00 1dy ... R y2 y1 KB 00 N+1dy 3775 and WN (t) = [wN 1 (t); ; wN N+1(t)]T contains the generalized Fourier coe cients. The approximate solution at time t is determined by solving the ODE system to that point and constructing the solution through the expansion (4.10). We reiterate that the coe cients ~ (y);g EI(y); g cDI(y) are spatially varying due to the nonhomogeneities from the patches and damage. 4.4 Finite Element Analysis A second approximation method which is widely used in structural applications is the nite element method. While this latter method shares a common mathematical background with 11 the previously described spline-based Galerkin technique, the manner of application is su ciently di erent to warrant discussion. We present the technique from two perspectives. The rst is formal and is intended to illustrate the construction of the resulting matrix system and highlight issues which must be addressed when employing the nite element method in damage detection applications. Once the method has been illustrated in this manner, we place it in a rigorous mathematical framework through a brief discussion of the Ciarlet triple and appropriate approximation spaces. The basic strategy in the nite element method is to consider the structure as an assemblage of individual elements as depicted in Figure 4. Associated with a typical element are nodal points at which we will de ne parameters which ultimately de ne the beam displacement. In this case, each element has two nodal points which are located at the element ends. Possible motions at the nodes are described as local degrees of freedom (local DOF). For the beam elements, displacements and rotations (slopes) are considered at each end of the element. This yields a total of 4 local DOF for the element. The displacements result from nodal forces whereas nodal moments yield the rotations. Note that by specifying both displacements and slopes, we will obtain a Hermite nite element as compared with Lagrangian elements in which only displacements are speci ed. c w(t,y) y y y 1 2 Figure 4. Finite element discretization of the cantilever beam. To uniquely determine the displacements and slopes at the element edges (4 local DOF), polynomials of at least degree 3 are required. These polynomials are typically de ned on a local interval K with local nodal values matched to obtain global displacements. To illustrate, we consider an element of length L located between `1 and `2; hence K = [0; L]. On this interval, the local (element) de ections are taken of the form we L(t; y) = a0(t) + a1(t)y + a2(t)y2 + a3(t)y3 : (4.6) With the notation fA(t)gT = [a0(t); a1(t); a2(t); a3(t)] and fz(y)g = [1; y; y2; y3]T , the approx12 imate solution and squared strain can be expressed as we L(t; y) = fA(t)gTfz(y)g @2we L @y2 (t; y)!2 = fA(t)gTfz00(y)gfz00(y)gTfA(t)g = fA(t)gT [D(y)]fA(t)g ; where [D(y)] fz00(y)gfz00(y)gT is a 4 4 matrix. Moreover if we de ne w1(t) = we L(t; 0) ; w2(t) = we L(t; L) 1(t) = @we L @y (t; 0) ; 2(t) = @we L @y (t; L) then enforcement of (4.11) yields fw(t)gi = 26664 w1(t) 1(t) w2(t) 2(t) 37775 = [B]fA(t)g ; [B] = 26664 1 0 0 0 0 1 0 0 1 L L2 L3 0 1 2L 3L2 37775 where fw(t)gi contains the nodal displacements and slopes. Equivalently, with C de ned as B 1, the original coe cients can be expressed as fA(t)g = [C]fw(t)gi. This then implies that we L can be represented as we L(t; y) = fz(y)gT [C]fw(t)gi : (4.7) To determine the motion, Hamilton's principle is again invoked. Demonstrating rst with the conservative case, the potential and kinetic energies for the element are expressed as U = 1 2 Z L 0 g EI(y) @2we L @y2 (t; y)!2 dy = fw(t)gTi [C]T 1 2 Z L 0 g EI(y)[D(y)]dy [C]fw(t)gi and T = f _ w(t)gTi [C]T 1 2 Z L 0 e (y)[F (y)]dy [C]f _ w(t)gi where [F (y)] fz(y)gfz(y)gT . Hamilton's principle then yields [M e]f  w(t)gi + [Ke]fw(t)gi = 0 (4.8) where [M e] = [C]T Z L 0 e (y)[F (y)]dy [C] [Ke] = [C]T Z L 0 g EI(y)[D(y)]dy [C] denote the local mass and sti ness matrices. 13 The vector equation (4.13) describes the dynamics of the unforced and undamped local element on [`1; `2]. To obtain global equations, the local dynamics are matched through enforcement of compatibility conditions regarding displacements and rotations (slopes). Speci cally, displacement and rotations at the element nodes are required to match those of neighboring elements. This is analogous to noting that moments and forces at the adjacent nodes are balanced. To illustrate, we consider a simpli ed case in which a homogeneous beam (hence constant parameters) is discretized into two identical elements. The local sti ness matrix in this case has the structure [Ke] = EI L3 26664 12 6L 12 6L 6L 4L2 6L 2L2 12 6L2 12 6L 6L 2L2 6L 4L2 37775 : Enforcement of displacement and rotational components yields the global sti ness matrix [K] = EI L3 2666666664 12 6L 12 6L 0 0 6L 4L2 6L 2L2 0 0 12 6L 24 0 12 6L 6L 2L2 0 8L2 6L 2L2 0 0 12 6L 12 6L 0 0 6L 2L2 6L 4L2 3777777775 : Note that in constructing the global system, the number of global DOF is the total number of local DOF minus the number of enforced compatibility conditions. The nonconservative damping and external e ects are incorporated via an extended Hamilton's principle as mentioned in Section 4.2. We emphasize that when considering the beam having surface-mounted piezoceramic patches and missing material at yc, the density, sti ness and damping parameters are not constant and quadrature rules must be employed when constructing the local mass, sti ness and damping matrices. Due to the piecewise nature of the parameters, the elements must be aligned with patch edges to preserve the accuracy of the method and care must be exercised when enforcing compatibility constraints. In summary, routines which admit spatial variability (including piecewise discontinuities) must be employed to obtain success when employing nite elements in damage detection schemes using piezoceramic actuators and sensors. 4.4.1 Mathematical Formulation of the Finite Element Method To place the method in the Sobolev framework outlined in Section 2, it is advantageous to consider a slightly more abstract de nition of a nite element and using this, construct approximating subspaces. Following the de nition of Ciarlet [14], a nite element is de ned as the triple (K;P;N ) where K denotes a local geometrical region, P is a nite dimensional space of functions on K, and N is the degrees of freedom [13, 14]. For this discussion, K = [0; L] [0; `] is an interval, P = P3 denotes the set of cubic polynomials and N = fN1; N2; N3; N4g where N1(v) = v(0); N2(v) = v0(0); N3(v) = v(L); N4(v) = v0(L) for all v 2 P. Note that N forms a basis for the dual space P . 14 For the nite element (K;P;N ), we let f 1; 2; 3; 4g denote a basis for P = P3 on the interval K. In terms of the previous discussion, this basis, termed the nodal basis, is given by fzgT [C] for K = [0; L] (see (4.12)). For L = 1, the four basis functions are plotted in Figure 5. Note that Ni( j) = ij for i; j = 1; ; 4. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y−axis Psi_1(y) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y−axis Psi_2(y) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y−axis Psi_3(y) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 y−axis Psi_4(y) Figure 5. Local Hermite basis functions 1(y); 2(y); 3(y); 4(y) de ned on the interval [0; L] = [0; 1]. The local interpolant is then de ned by the expansion IKv 4 X j=1Ni(v) i (see [13, 14] for the general formulation of the local interpolant). Note that in the notation of the previous discussion, the nodal values are given by fwig = [w1; 1; w2; 2]T and the local interpolant is represented by we L = fzgT [C]fwgi (again see (4.12)). With local interpolants thus de ned, it is necessary to de ne a global interpolant over the full domain. To this end, let T denote a subdivision of [0; `] into a set of subintervals Kj with [0; `] = SKj . The global interpolant IT is then de ned as IT gjKj = IKjg 15 for g 2 C1[0; `] and Kj 2 T . By matching nodal values (hence displacements and slopes)between elements, IT g 2 C1 for all g 2 C1[0; `] and the space V N = Vh C1[0; `] de ned byVh = fg : gjKj 2 P3(Kj) for all Kj 2 T ; and g; g0 continuous at the nodesgis termed a C1 nite element space (the notation Vh is common in the nite elements liter-ature and we retain it here to denote the approximating nite dimensional space). It is thisspace Vh which is considered in the approximation framework of Section 2. The reader isreminded that the regularity of the nite element space depends on the dimension and whilethe Hermite elements described here yield C1 interpolants on the interval, they will yield onlyC0 interpolants when triangulating in lR2 (see [13] and [15, pages 80,81]).The nite element approach di ers from the previously discussed cubic spline Galerkinmethod in the local-to-global nature of the approximation. Through an a ne mapping ofbasis functions on a reference element [0; L] to [`1; `2], local interpolants are constructed. Likethe cubic spline approximates, these local interpolants satisfy a smoothness criteria dictatedby the choice of basis. Global interpolants are then obtained though enforcement of nodalcompatibility criteria. This yields a C1 method satisfying mixed smoothness and interpolatorycriteria. The requirement of matching nodal displacements and slopes is unnecessary with thecubic spline basis since all smoothness properties are directly constructed in the basis. Thisyields a C2 approximate solution which satis es purely interpolatory constraints.5 ConclusionThis brief review has centered on a nondestructive evaluation (NDE) technology based uponPDE models of the physical system of interest. By formulating the method in an abstractoperator setting, one obtains well-posedness and a convergence theory which is su cientlygeneral to include NDE applications ranging from thermal imaging to vibration analysis withsmart material sensors and actuators.A key component in the implementation of the method is the choice of approximationscheme when discretizing the PDE model to obtain a matrix system. Due to the local nature ofactuator inputs and damage, approximation methods employing locally de ned basis functionsappear to yield superior sensitivity when solving the inverse problem to estimate physicalparameters. As demonstrated in investigations referenced in [2, 10], methods based uponglobally de ned modal bases and frequency shifts lack the sensitivity to adequately locateand characterize damage in general applications (the sensitivity often lies below the thresholdof measurement error). On the other hand, cubic spline based Galerkin methods have beensuccessfully employed to characterize damage representing less than 0:1% of the undamagedstructural state. Due to the issue of sensitivity with respect to local structural variations, it isalso hypothesized that Galerkin expansions employing globally de ned Legendre or Chebyshevbases will prove less successful in NDE applications than those employing spline or niteelement bases.The success of either the spline-based Galerkin methods or nite element approaches iscontingent upon the construction of approximate solutions which admit variable physicalparameters (including piecewise discontinuities). To retain accuracy with either approach,gridpoints must be aligned with the discontinuities (e.g., patch edges). This can be a limitation16 for codes employing automatic mesh-independent error analysis. The same issue must be facedif parameter discontinuities occur due to the nature of the damage. The determination of anappropriate mesh is more di cult in this latter case, however, since the region of discontinuityis unknown and is to be determined through parameter estimation.The models and numerical methods must also be exible with regards to boundary condi-tions. While this issue may at rst appear unrelated to the ultimate goal of damage detection,it is in fact crucial for attaining su cient accuracy to detect material nonhomogeneities. Bothnumerical investigations [9] and experiments [1] have illustrated the necessity of modelinginexact boundary conditions for structures with boundary energy loss as compared with theoption of compensating through modi cation of structural parameters (e.g., increased anddecreased E will decrease some frequencies in a manner similar to loosened boundary clamps).As illustrated in [1, 9], however, this latter technique does not provide accurate model tsover a signi cant frequency range nor does it lead to consistency among estimated parame-ters. From an NDE perspective, it is awed since the modi cation of physical parameters toaccount for boundary inaccuracies will tend to reduce the accuracy for determining variationsdue to damage. Hence numerical methods (either spline or nite element) must be su cientlyexible to approximate boundary variations which are incorporated in the PDE model.Acknowledgements This research was supported in part by the Air Force O ceof Scienti c Research under grants AFOSR-F49620-93-1-0198 and AFOSR-F49620-95-1-0236,by NASA under grant NAG-1-1600 and by the U.S. Department of Education through aGAANN Fellowship to Y.Z.References[1] H.T. Banks, D.E. Brown, R.J. Silcox, R.C. Smith and Y. Wang, Modeling and estimationof boundary parameters for imperfectly clamped structures, to appear.[2] H.T. Banks, D.J. Inman, D.J. Leo and Y. Wang, An experimentally validated damagedetection theory in smart structures, CRSC Technical Report CRSC-TR95-7, January,1995; Journal of Sound and Vibration, to appear.[3] H.T. Banks and K. Ito, A uni ed framework for approximation in inverse problems fordistributed parameter systems, in Control: Theory and Advanced Technology, 4(1), 1988,pp. 73-90.[4] H.T. Banks and F. Kojima, Approximation techniques for domain identi cation in twodimensional parabolic systems under boundary observations, Proc. 26th IEEE Conf. onDec. and Control, December, 1987, Los Angeles, 1411-1416.[5] H.T. Banks and F. Kojima, Boundary shape identi cation problems in two dimensionaldomains related to thermal testing of materials, Quart. Applied Math., 47, 1989, pp.273-293.[6] H.T. Banks, F. Kojima and W.P. Winfree, Boundary estimation problems arising inthermal tomography, Inverse Problems, 6, 1990, pp. 897-921.17 [7] H.T. Banks and K. Kunisch, Estimation Techniques for Distributed Parameter Systems,Birkhauser Boston, 1989, 320 pp.[8] H.T. Banks, R.C. Smith and Y. Wang, The modeling of piezoceramic patch interactionswith shells, plates and beams, Quart. Applied Math., 53, 1995, pp. 353-387.[9] H.T. Banks, R.C. Smith and Y. Wang, Parameter estimation for an imperfectly clampedplate { numerical examples, Proceedings of the 1995 Design Engineering Technical Con-ferences, Volume 3, Part C, Boston, MA, September 17-20, 1995, 963-972.[10] H.T. Banks, R.C. Smith and Y. Wang, Smart Material Structures: Modeling, Estimationand Control, CRSC Lecture Notes CRSC-LN96-1, March, 1996; Masson/J. Wiley, toappear.[11] H.T. Banks and Y. Wang, Damage detection and characterization in smart materialstructures, CRSC Technical Report CRSC-TR93-17, November 1993; in Control andEstimation of Distributed Parameter Systems; Nonlinear Phenomena, Birkhauser ISNM,118, 1994, pp. 21-43.[12] H.T. Banks, Y. Wang and D.J. Inman, Bending and shear damping in beams: frequencydomain estimation techniques, ASME J. Vibration and Acoustics, 116, 1994, pp. 188-197.[13] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods,Springer-Verlag, New York, 1994.[14] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, NewYork, 1978.[15] C. Johnson, Numerical Solution of Partial Di erential Equations by the Finite ElementMethod, Cambridge University Press, New York, 1987.[16] P.M. Prenter, Splines and Variational Methods, Wiley, New York, 1975.[17] R. Weinstock, Calculus of Variations with Applications to Physics and Engineering,Dover, New York, 1974.[18] J. Wloka, Partial Di erential Equations, Cambridge University Press, Cambridge, 1987.18